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The time-dependent numerical renormalization-group approach (TD-NRG), originally devised for 
tracking the real-time dynamics of quantum-impurity systems following a single quantum quench, 
is extended to multiple switching events. This generalization of the TD-NRG encompasses the 
possibility of periodic switching, allowing for coherent control of strongly correlated systems by an 
external time-dependent field. To this end, we have embedded the TD-NRG in a hybrid framework 
that combines the outstanding capabilities of the numerical renormalization group to systematically 
construct the effective low-energy Hamiltonian of the system with the prowess of complementary 
approaches for calculating the real-time dynamics derived from this Hamiltonian. We demonstrate 
the power of our approach by hybridizing the TD-NRG with the Chebyshev expansion technique 
in order to investigate periodic switching in the interacting resonant-level model. Although the 
interacting model shares the same low-energy fixed point as its noninteracting counterpart, we sur- 
prisingly find the gradual emergence of damped oscillations as the interaction strength is increased. 
Focusing on a single quantum quench and using a strong-coupling analysis, we reveal the origin of 
these interaction-induced oscillations and provide an analytical estimate for their frequency. The 
latter agrees well with the numerical results. 

PACS numbers: 03.65.Yz, 73.21. La, 73.63.Kv, 76.20.+q 



I. INTRODUCTION 

The quantitative description of real-time dynamics in 
strongly correlated systems is one of the outstanding 
challenges of contemporary condensed-matter physics, 
with relevance to varied systems ranging from cold 
atomai^ and dissipative systems^ to quantum-dot de- 
vices^ and biological donor-acceptor molecules^ Along- 
side fundamental question concerning the underlying 
time scales and the long-time behavior, there are many 
technological issues that require careful investigation. 
For example, the decoherence and relaxation of spins 
appears to be the major obstacle for the realization 
of quantum-computing algorithms in real systems^ An- 
other key issue is the understanding of coherent control 
and the switching characteristics of nanodevices such as 
single-electron transistors^ These and related topics re- 
quire the development and application of suitable many- 
body techniques^ 

Over the years, the Kadanoff-Baym 10 and Keldyshii 
techniques have proven to be accurate tools for describ- 
ing the real-time dynamics of weakly correlated systems 
such as light-matter interaction in semiconductors^ 2 - and 
the decoherence and relaxation of an impurity spin well 
above the Kondo temperature^ Geared toward pertur- 
bation theory, these techniques generally fail upon the 
development of strong correlations, when nonperturba- 
tive approaches are in order. A case in point are quan- 
tum dots tuned to the Kondo regime,— where traditional 
diagrammatic-based approximations are unsuitable to 
describe the nonequilibrium stated The difficulty lies 



in the fact that strongly correlated systems change their 
nature as a function of certain control parameters such 
as the temperature or the coupling constants, an aspect 
well captured by renormalization-group approaches ! 16 > 17 
The precise status of a voltage bias as yet another con- 
trol parameter in interacting nanostructures is still under 
debate. 

Recent years have witnessed an impressive advance- 
ment of numerical techniques aimed at tracking the 
real-time dynamics of strongly correlated systems, pri- 
marily with the development of the time-dependent 
density-matrix renormalization group (TD-DMRG) »i2r— 
Yet while the adaptive TD-DMRG works remarkably well 
on time scales of order the reciprocal bandwidth, it is 
presently unsuited for tackling longer time scales due to 
an accumulated error that grows first linearly and then 
exponentially with the time elapsed. Alternative formu- 
lations 2 - of the TD-DMRG circumvent the accumulated 
error, but are simply too demanding to advance to long 
times. Recent adaptations of continuous-time Monte 
Carlo techniques to real-time dynamics 2 ^— are free of 
finite-size effects, but are confined to short time scales 
due to an inherent sign problem. The Chebyshev expan- 
sion technique^ 7 - developed by Tal Ezer and Kosloff ) 28 i 29 
offers yet another extremely powerful approach for track- 
ing the time evolution of finite-size systems. However, 
it too is quite limited in accessing long time scales in 
the presence of interactions due to the exceedingly large 
Hilbert space that must be retained. A complementary 
approach is provided by the time-dependent numerical 
renormalization group (TD-NRG) The TD-NRG 
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can successfully bridge over vastly different time scales, 
but is far more restrictive in the systems and perturba- 
tions to which it can be applied. A composite approach 
that combines the complementary traits of the different 
techniques mentioned above is highly desirable. 

In this paper, we devise such a hybrid approach that 
combines the outstanding capabilities of the numeri- 
cal renormalization group (NRG) to systematically con- 
struct the effective low-energy Hamiltonian of the sys- 
te m 16 ' 33 with the prowess of complementary approaches 
for calculating the real-time dynamics derived from this 
Hamiltonian. Typically strongly correlated systems pos- 
sess multiple energy scales that markedly differ in mag- 
nitude, hence their dynamics is governed by vastly dif- 
ferent time scales. This spread of time scales, which may 
differ by many orders of magnitude, poses an enormous 
obstacle for most computational approaches. The TD- 
NRG is quite unique in this respect as it can efficiently 
bridge between the different time scales. Our hybrid ap- 
proach presented below provides a flexible platform for 
combining the TD-NRG with one's method of choice for 
treating the effective low-energy Hamiltonian. Possible 
choices for the complementary method could be exact 
diagonalization, the TD-DMRG, and possibly also the 
time-dependent noncrossing approximation! 34 ' 35 Here we 
shall demonstrate the applicability of our approach by 
combining the TD-NRG with the Chebyshev expansion 
technique (CET), providing thereby an important proof 
of principle. 

The basic philosophy underlying the hybrid-NRG is 
to first exhaust the TD-NRG in order to decompose the 
time-dependent wave function into distinct components, 
each associated with a separate time scale and evolving 
according to its own reduced Hamiltonian acting on a 
suitable subspace of the full Fock space. Each of the in- 
dividual components with its associated Hamiltonian can 
then be treated with improved accuracy using, e.g., the 
TD-DMRG or CET. In this manner one can exploit the 
successive reduction in energy scales in order to boost the 
TD-DMRG and CET to long time scales that otherwise 
would be inaccessible to either of these methods. Con- 
comitantly, the accuracy and flexibility of the TD-NRG 
are greatly enhanced, as we demonstrate by extending 
the approach to the physically relevant case of repeated 
switchings. The hybrid platform further offers an ap- 
pealing way to reduce discretization errors inherent to 
the Wilson chain by converting to hybrid chains. 

The idea to use the NRG level flow to construct effec- 
tive Hamiltonians is, of course, not new. Dating back to 
the original work of Wilson^— this framework has been 
significantly advanced by Hewsor*^ 3 . who used it as a 
starting point for devising a renormalized perturbation 
theory. Hewson's approach requires, however, analytical 
knowledge of the form of the low-energy Hamiltonian and 
its associated quantum field theory. So far it has been 
applied mainly to the single-impurity Anderson model, 
where it was used, among other things, to calculate the 
steady-state curren t 33 ' 36 in the limit of a small bias volt- 



age. 

Our approach is far more general as it makes no as- 
sumption on the analytical form of the effective Hamil- 
tonian. Rather, it is solemnly based on Wilson's original 
concept^ 6 - that the NRG level flow contains all accessible 
information, and in particular can accurately describe 
the crossover region between two distinct fixed points (a 
regime which generally lies outside the reach of pertur- 
bative methods). Our framework exclusively uses the 
sequence of diagonalized NRG Hamiltonians ^ circum- 
venting thereby any prejudice on the form of the effective 
Hamiltonian. As a result our method is model indepen- 
dent, relying solely on the NRG approach itself. 

In this paper, we extend the original TD-NRG algo- 
rithm from a single quantum quench to multiple switch- 
ings, which requires an additional approximation beyond 
the one underlying the conventional TD-NRG. As we 
demonstrate by explicit calculations, the quality of the 
approximation can be systematically improved by en- 
larging the subspace treated using the complementary 
method. 



A. Preliminaries 

In the TD-NRG, the continuous bath is represented 
by a discretized Wilson chain^ characterized by tight- 
binding hopping matrix elements that decay exponen- 
tially along the chain. This separation of energy scales 
enables access to exponentially long time scale a 30 ' 31 that 
otherwise would be inaccessible using ordinary tight- 
binding chains. In analyzing the accuracy of the TD- 
NRG it is important to distinguish between two sources 
of error: one extrinsic due to the discretized representa- 
tion of the continuous bath in terms of a Wilson chain, 
and the other intrinsic due to the TD-NRG algorithm for 
tracking the real-time dynamics on the Wilson chain. 

As already noted in Ref. Hi], the latter source of error is 
remarkably small. We illustrate this point in Fig.[lja) for 
the noninteracting resonant-level model (RLM) , describ- 
ing a single fermionic level coupled by hybridization to a 
conduction band (see Sec. [V] for an explicit definition of 
the model). Abruptly shifting the energy of the level and 
tracking the time evolution of the level occupancy rid(t), 
we compare the TD-NRG to an exact analytical solu- 
tion for a continuous bath [given by Eq. (50) of Ref. l3l| . 
as well as to an exact numerical solution on the Wilson 
chain using exact diagonalization of the single-particle 
eigenmodes. Only minuscule deviations are found be- 
tween the TD-NRG and the exact solution on the Wilson 
chain, both of which significantly depart at some point 
from the continuum-limit result. As analyzed in Ref. l3ll . 
one can decrease the deviations from the continuum limit 
and delay them to a later time by reducing the Wilson 
discretization parameter A. At the same time, the devi- 
ations are hardly affected by prolonging the chain length 
N. 

This leads to two important conclusions: (i) The main 
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FIG. 1. (Color online) (a) Real-time dynamics of the impurity 
charge n c i{i) in the resonant-level model following a sudden 
quench of the level energy from Ed = — 2Fo to Ed = 2Fo at 
time t = 0. The red solid line shows the exact time evolution 
on the Wilson chain, obtained by exact diagonalization of the 
single-particle eigenmodes. The dashed blue line depicts the 
TD-NRG result, while the solid black line shows the exact 
analytical solution for a continuous band [given by Eq. (50) 
of Ref. . (b) A two-dimensional contour plot of the exact 
time-dependent occupancies m (t) of the first 26 sites along the 
Wilson chain, labeled i = 0, . . . , 25. The red arrows indicate 
instances in time when reflected currents reach the impurity 
site. The very same times are marked by the black arrows in 
panel (a). Parameters: A = 1.8, Yq/D = 10 -2 , the number of 
states kept in the TD-NRG is N s — 800, and the chain length 
is N = 40. 



source of error in the TD-NRG is extrinsic rather than 
intrinsic^ (ii) The exponentially decaying tight-binding 
matrix elements, which lie at the heart of the NRG,— are 
also the limiting factor for reproducing the continuum- 
limit result for the quench dynamics. 

To understand the source of the deviations from the 
continuum-limit result, we note that globally conserved 
quantities such as the charge or the spin of the system 
are locally connected by the continuity equation to as- 
sociated charge and spin currents. The exponentially 
decreasing tight-binding matrix elements along the Wil- 
son chain significantly slow down the propagation of such 
currents^ generating internal reflections at the sites of 
the one-dimensional chain. This is depicted in Fig. [TJb) , 
where we plot the exact time-dependent occupancies of 
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FIG. 2. (Color online) The exact single-particle energy levels 
of the RFM on a Wilson chain with A = 1.8, V /D = 10" 2 , 
and Ed = 0. The chain length is N = 40. Panel (a) depicts 
positive energies on a logarithmic scale. Panel (b) shows the 
full spectrum on a linear scale. 



the first 26 sites along the chain in response to a sud- 
den quench of the impurity level. The two-dimensional 
contour plot clearly reveals reflections of the charge cur- 
rent at certain positions and certain characteristic times. 
Once the reflected currents reach the impurity site its 
level occupancy starts to deviate significantly from the 
exact continuum-limit result. Indeed, the red arrows in 
Fig. [IJb) indicate instances in time when reflected charge 
wave fronts reach the impurity site. At these very same 
times the occupancy on the impurity level develops new 
structures [marked by the black arrows in Fig. Ufa)] that 
are absent in the continuum limit. Upon decreasing the 
Wilson discretization parameter A the magnitude of the 
reflected currents is suppressed and the reflection points 
are pushed deeper down the chain, however the effect is 
never fully eliminated as long as A > 1. 

An alternative perspective on the effect of the Wilson 
discretization procedure is provided by examining the ex- 
act single-particle energy levels of the RLM. For an or- 
dinary tight-binding chain of length N with a constant 
hopping matrix element £, the single-particle energy lev- 
els are roughly uniformly distributed in the energy range 
[-D, D], where D = 2£ is the conduction-electron band- 
width. As the chain length is increased the energy levels 
become more densely distributed until a continuous spec- 
trum is recovered for N — > 00. A different picture applies 
to the Wilson chain. As depicted in Fig. the single- 
particle energy levels are uniformly distributed on a loga- 
rithmic scale, resulting in a sparse distribution of levels at 
higher energies. Enumerating the positive single-particle 
energy levels from high to low, these scale as e„ cx A~ n , 
in accordance with Wilson's logarithmic discretization of 
the conduction band^ By prolonging the chain length N 
one increases the total number of levels, yet the distribu- 
tion of high-energy levels retains its sparse form even as 
A^ — > 00. A continuous band is recovered only upon im- 
plementing the combined limit A — > 1 + , A" — > 00, which 
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illustrates the limitation of working with a fixed A > l. 38 
These argumentations clearly point to an intrinsic 
tradeoff within the TD-NRG, as the very same logarith- 
mic discretization that enables access to exponentially 
long time scales also prevents from fully recovering the 
continuum limit. The question then arises whether one 
can somehow reconcile these two apparently contradict- 
ing properties, which is one of the main goals of this 
work. 



B. Plan of the Paper 

Briefly stated, the aim of this paper is a two-fold ex- 
tension of the TD-NRG. 

The first goal is to devise a flexible framework for hy- 
bridizing the TD-NRG with complementary methods of 
calculating real-time dynamics that do not rely on the 
special structure of the Wilson chain. By liberating our- 
selves from working on the Wilson chain at all time scales, 
the hybrid method should enable a systematic improve- 
ment of both finite-size and discretization errors. 

The second objective is to extend the method from a 
single quantum quench to repeated switchings between 
two distinct Hamiltonians H a and W b . Besides being of 
considerable interest on its own right, this protocol can 
be viewed as a first vital step toward treating a general 
time-dependent Hamiltonian H(t). Indeed, discretizing 
the time axis and replacing H(t) within each time interval 
i with the constant form "Hi, the full time evolution can 
be interpreted as a sequence of quenches. 

Despite its conceptual simplicity this strategy has not 
been pursued thus far since the initial nonequilibrium 
density operator pi at the beginning of each time inter- 
val is neither explicitly known nor of the analytical form 
required by the original TD-NRG formulation^ 

To this end, we derive the hybrid-NRG in Sec. [TIJ Since 
the formulation relies heavily on the TD-NRG, we shall 
commence in Sec. Ill Al with a brief review of the TD- 
NRG following Ref. 31. This exposition is essential not 
only for keeping the paper self-contained, but mainly for 
introducing the notations and tools that will be used 
throughout our construction of the hybrid-NRG. After 
setting the stage in Sec. Ill Al we present the first major 
conceptual result of this paper in Sec. Ill Bl the hybrid- 
NRG approach. The key idea is to partition the Wil- 
son chain into two parts: the high-energy part is treated 
with the TD-NRG while the low-energy part is feed into 
the complementary method of choice. The standard TD- 
NRG approach can be embedded into this more general 
hybrid framework by shifting the partition to the end of 
the Wilson chain. In order to make contact with the com- 
plementary approach used for hybridizing with the TD- 
NRG (e. g., the CET), we examine the interface with the 
TD-NRG from a wave-function perspective in Sec. lHB^I 

Using the central results of Sec. UTJ we extend the hy- 
brid approach to periodic switching between two Hamil- 
tonians in Sec. Mil making the simulation of coherent con- 



trol by an external field accessible to the hybrid method. 
This section covers the second major conceptual result 
of this paper. Since the exact evaluation of the reduced 
density matrix of an arbitrary density operator has re- 
mained a computational challenge, we propose a set of 
approximations which neglect some of the high-energy 
contributions to the real-time dynamics. These approxi- 
mations are systematic and are carefully analyzed. The 
complementary method used to supplement the TD-NRG 
in this paper, the Chebyshev expansion technique , 28 i 29 is 
reviewed in Sec. IIVI 

Since the deviations of finite-size nonequilibrium dy- 
namics from the continuum limit are linked to the bath 
discretization, our hybrid framework targets the libera- 
tion of numerical simulations from the particular form of 
the Wilson chain without losing access to exponentially 
long time scales. In Sec.fV] we demonstrate the potential 
of our hybrid framework by investigating the influence of 
different hybrid chain types on the discretization errors 
encountered in the local dynamics. In Sec. lVIl we present 
our results for periodic switching. Since the TD-NRG is 
embedded in our hybrid framework, we first discuss in 
Sec. IVI A II the limitations of a simple extension of the 
TD-NRG to periodic switching. The exact solution of 
the finite-size RLM subject to a periodic drive is used to 
benchmark both the hybrid NRG-CET and the periodic 
TD-NRG. 

The interacting resonant-level model 3 ^— (IRLM) 
serves as a first nontrivial test of our hybrid approach. 
It includes an additional local capacitive coupling U be- 
tween the charge on the impurity level and the fermionic 
band, preventing an exact solution of its dynamics. The 
low-energy fixed point of the IRLM is identical to that 
of the noninteracting RLM, featuring a renormalized, 
{/-dependent hybridization strength. By comparing the 
real-time dynamics of both models with identical renor- 
malized hybridization strengths, we show in Sec. IVI Bl 
that the local charge dynamics significantly deviates with 
increasing U from the noninteracting case. Interaction- 
induced oscillations are found in the IRLM, whose char- 
acteristic frequency depends only on the renormalized 
hybridization and not directly on U. Those oscillations 
are completely absent in the noninteracting RLM, even 
though both models share the same low-energy fixed 
point. Using a strong-coupling analysis, we provide a 
simple physical picture for this surprising effect. Finally, 
we conclude with a discussion and outlook in Sec. IVIII 



II. THE HYBRID-NRG 

A. The time-dependent NRG 

The TD-NRG has been designed to track the real- 
time dynamics of quantum-impurity systems following an 
abrupt quantum quench. The perturbations under con- 
sideration are implicitly assumed to be of local character, 
i.e., perturbations that are applied either to the impurity 
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FIG. 3. The full Wilson chain of length TV is divided into 
a subchain of length m and the "environment" R m ,N- The 
Hamiltonian H m can be viewed either as acting only on the 
subchain of length m, or as acting on the full chain of length 
TV, but with the hopping matrix elements t m ,--- , ijv-i all 
set to zero. The former picture is the traditional one. In the 
TD-NRG one adopts the latter point of view. 

itself or to its close vicinity. 

The Hamiltonian of a quantum-impurity system has 
the generic structure 

H = %bath + %imp + %mix, (1) 

where 'Hbath models the continuous bath, Hi mp represents 
the decoupled impurity, and H m i x describes the coupling 
between the two subsystems. The entire system is char- 
acterized at time t = by the density operator 

e -f3W 

k = Trace {e^}' (2) 

when a static perturbation A% is suddenly switched on: 
U(t > 0) = W+AU = H f . The density operator evolves 
thereafter in time according to 

p(t>0) = e-^ f p e mf . (3) 

Our objective is to use the NRG to compute the time- 
dependent expectation value 0(t) of a general local op- 
erator O. As shown in Ref. and|3l|, the result can be 
written in the form 

N trun 

= E E j t{E ™- ET) 0™pl e *(m), (4) 

m r,s 

where E™ and E™ are the dimension-full NRG eigenen- 
ergies of the perturbed Hamiltonian at iteration to < N, 
0™ s is the matrix representation of O at that iteration, 
and p r s c ^(m) is the reduced density matrix defined in 
Eq. (|10[) below. The restricted sum over r and s requires 
that at least one of these states is discarded at iteration 
to. The NRG chain length TV implicitly defines the tem- 
perature entering Eq. @: T/v oc A _JV / 2 , where A > 1 is 
the Wilson discretization parameter. 

The derivation of Eq. (|4]) relies on two key ingredients: 
(i) The identification of a complete basis set of approxi- 
mate NRG eigenstates for the many-body Fock space Fn 
of the Wilson chain; (ii) Expectation values are obtained 
by explicitly tracing over this complete basis set using a 
suitable resummation procedure. Below we review these 
two key components following the notations and presen- 
tation of Ref. l3~ll 



1. Complete basis set 

The NRG targets an iterative solution of a quantum 
impurity coupled to a finite Wilson chain with TV chain 
links. 44 Similar to the initial sweep in the finite-size 
DMRG, one can view the NRG procedure as a set of op- 
erations, where at first all hopping matrix elements are 
set to zero along the TV-site chain, and at each successive 
step another hopping matrix element is switched on. The 
full Hamiltonian "Hat is recovered once all hopping matrix 
elements have been switched on. The entire sequence of 
Hamiltonians H m with m < TV act in this picture on the 
same Fock space Fn of the TV-site chain, hence each NRG 
eigenenergy of % m has an extra degeneracy of S- N ~ m \ 
where d is the number of distinct configurations at each 
site along the chain. The extra degeneracy stems from 
the TV — to "environment" sites at the end of the chain, 
denoted by R m ,N in Fig- El which remain decoupled from 
the impurity at iteration m. 

When acting on the TO-site chain, we label the NRG 
eigenstates and eigenenergies of H m by {|r;m)} and 
E™, respectively. Consider now the action of H m on 
the full TV-site chain. Enumerating the different con- 
figurations of site i by {ai}i=i....,d, each of the tensor- 
product states |r;m) (g) |a m +i, . . . , a/v) with arbitrary 
a m+ i, . . . , ajv is a degenerate NRG eigenstate of Hm with 
energy £™ . To label these states we introduce the short- 
hand notation |r, e;m), where the "environment" vari- 
able e = {a m +i, . . . , a^v} encodes the TV — m site labels 
a m+ i, . . . , otjv, and the index to is used to record where 
the chain is partitioned into a "subsystem" and an "en- 
vironment" (see Fig. [3]). 

In order to retain a manageable number of states, the 
high-energy eigenstates are discarded after each iteration, 
which is fully justified in equilibrium by the hierarchy of 
energy scales along the Wilson chain and the Boltzman- 
nian form the equilibrium density operator. Regarding 
all states of the final iteration as discarded, it has been 
show in Refs. and [3l| that the collection of all states 
discarded in the course of the NRG iterations form a com- 
plete basis set of approximate NRG eigenstates for the 
full TV-site chain. 

To understand this important point, consider the first 
iteration TO m i n at which states are discarded. In or- 
der to keep track of the complete basis set of the TV- 
site chain, the eigenstates |r, e; TO m i n ) can be formally 
divided into two distinct subsets: the discarded high- 
energy states {|/, e; m m i n ) t ji s } and the kept low-energy 
states {|fc, e; TO m i n )kp}. Obviously, the sum of the two 
subsets form a complete basis set of the full chain. To 
simplify the notations we shall omit hereafter the sub- 
scripts dis and k p , and will use in exchange the indices I 
and k to label the discarded and kept states, respectively. 
At the next NRG iteration only the kept states are used 
to construct the NRG eigenstates of %m min +i within the 
truncated subspace spanned by {\k, e; TO m j n )}. The re- 
sulting NRG eigenstates can again be divided into two 
subsets of discarded and kept states which, when com- 
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bined with the discarded eigenstates of iteration m m i n , 
form a complete basis set of the Fock space J-jy of the 
full JV-site chain. Repeating this procedure at all sub- 
sequent iterations, one continues to maintain a complete 
basis set of Tn up to the final NRG iteration m = N. 
In this manner we arrive at the following completeness 
relation 

N N 

E ^2\l,e;m)(l,e;m\= ^ P m = 1, (5) 

m— m m i n /,e ra— m m i n 

where the summation over to starts from the first iter- 
ation m m i n at which a basis-set reduction is imposed. 
Here the summation indices I and e implicitly depend 
on to, and the projector onto the subspace discarded at 
iteration m m i n < m < N is defined as 

P m =}^ \l, e; m){l, e; m\ . (6) 

l,e 

The complementary projector 1+ onto the subspace 
retained at iteration to (to < N) is given in turn by 

1«> = E \ k , e ; m )(k,e;m\ , (7) 

fc,e 

which can be recast in the form 

N 

m'— m+l 

This latter equality reflects the fact that all states re- 
tained at iteration m are necessarily discarded at some 
later iteration m' . In particular, since all states of the 
final iteration N are regarded discarded, then 1^ is iden- 
tically zero while l^r-i coincides with Pjy. Combined 
with Eq. ([5]), the completeness relation of Eq. (|S|) can 
further be partitioned into 

M 

1= £ Pm + itt (9) 

m=m m i n 

with arbitrary m m i n < M < N . This useful identity will 
be repeatedly used in constructing the hybrid- NRG. 



p r s c ^(m), defined as 

e 

Here the states \r,e;m) and \s,e;m) correspond to the 
Hamiltonian rU ', and the summation runs over the en- 
vironment degrees of freedom e. The only approxima- 
tion entering Eq. ^ is the standard NRG approximation 
Hn\t, e;m) H m \r,e;m) = E™\r, e;m), which enables 
us to write 

(s, e; m\p(t)\r, e; m) = e ltl ^ Er ~ E ° \s , e; m\po\r, e; m) . 

(11) 

Apart from this sole point, Eq. ^ constitutes an exact 
evaluation of 0(t) on the discretized iV-site chain. 

Practical calculations hinge on the ability to accurately 
compute the reduced density matrices of Eq. (fTUj) . For a 
general po this can be a daunting task. However, in the 
case of interest where po has the standard Boltzmann 
form of Eq. (J2J), the summation over e can be carried 
out exactly. Hence p T s e £ (m) can be evaluated at the same 
level of accuracy as the equilibrium density operator po- 
Technically this goal is achieved by implementing two in- 
dependent NRG runs, one for the initial Hamiltonian W 
in order to construct po using to Eq. and another 
for the full Hamiltonian . The reduced density matrix 
pl C r{m) is first evaluated with respect to the eigenstates 
of the initial Hamiltonian, and then rotate d 30 ' 31 to the 
eigenstates of the full Hamiltonian using the overlap ma- 
trices 

(q l ;m\r; m) = S qur (m). (12) 

Here |r; m) denotes an NRG eigenstate of the full Hamil- 
tonian at iteration to, and \qi]m) is an NRG eigenstate 
of the initial Hamiltonian at the same iteration. The 
method of calculating the overlap matrices S qi ^ r {m) is 
detailed in Ref. [3l|. 



B. Derivation of the Hybrid-NRG 



2. Reduced density matrix and the TD-NRG algorithm 

Following a sudden quench, the time evolution of the 
system is governed by the perturbed Hamiltonian W 
while the initial condition is encoded in the initial den- 
sity matrix pq. For quantum- impurity systems, all rele- 
vant information on the initial condition is contained in 
Eq. (H|) in the form of the reduced density matrice o 30 ' 31 



The original TD-NRG approach, summarized above, 
tracks the quench dynamics of a quantum-impurity sys- 
tem in terms of the phase factors e lt ^ Er ~ Es ' and the 
reduced density matrices p r s °^(m) assigned to each NRG 
iteration to. Although quite elegant and useful, it is less 
transparent how to incorporate ideas from methods such 
as the TD-DMRG or CET, as these deal with wave func- 
tions directly. To develop a convenient and flexible in- 
terface between the TD-NRG and these vastly different 
approaches we reformulate the former approach from a 
wave-function perspective. 
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1. Wave-function formulation 



Eq. (08]) as 



Let us commence with accurately stating the problem 
from a wave-function perspective. We are interested in 
tracking the time evolution of some initial state l^o) un- 
der the dynamics defined by the Hamiltonian % acting 
on a finite Wilson chain of length N. We shall not con- 
cern ourselves at this stage with how the initial state \4>o) 
is generated, but will elaborate on this important point 
later on. 



Formally our task boils down to computing 



mi 



(13) 



Application of the completeness relation of Eq. (J9j> to the 
state \ijj{t)) leads to its partitioning according to 



A I 



where 



and 



\^ m {t)) = p m m)) 



\xM(t)) = itt\m) 



(15) 



(16) 



are the projections of \ip(t)} onto the subspaces defined by 
P m and l~l l7 respectively. Equation (fT4|) simply converts 
the general state \if>(t)) into a concrete representation in 
terms of our complete basis set. 



E 



M 



\m'— m+1 

+ ( E Pm ' + *-it) Ap " 

\m' — m+1 / 



and noting that 



AI 



1+ - V P , 4- 1 + 

m'— m+1 

the operator A is recast in the exact form 

M 



A 



E M 



m) + A 



E \^(t)) + \XM(t)) , (14) where 



A{m) — P m AP m + Y^APjn + P m Al r 



and 



A — i + 41 + 



(19) 



(20) 



(21) 



(22) 



(23) 



Here the index M can take any value in the range m m j n < 
M < N. Explicitly, the operator A(m) has the formal 
representation 

trun 

M m ) = E E \r,e; rn)(r,e;m\A\s,e'; rn)(s,e';m\ , 

r,s e,e' 

(24) 

where the restricted sum J^/" 11 implies, as before, that 
at least one of the states r and s is discarded at iteration 
rn . 



2. Evaluation of expectation values 

Given the time-evolved wave function of Eq. (|14|) . we 
proceed to compute time-dependent averages of physical 
observables: 

A{t) = (m\Mm) ■ (17) 

To this end, we use the completeness relation of Eq. (|9]) 
to decompose any arbitrary operator A into 

M M 

A = ^ ] P m AP m i + ^ ^ ^i^ m Ai^[j + i^j74_P m |' 

m,m' m 

+ nA+M , (is) 

where the summations over m and m! start from m m i n . 
Writing the first two terms on the right-hand side of 



As in the original TD-NRG, we focus hereafter on lo- 
cal operators A that act solely on degrees of freedom that 
reside either on the impurity itself or on the first m m j n 
sites along the Wilson chain^L For any such local opera- 
tor, the matrix elements in Eq. (I24[) are diagonal in and 
independent of the environment variables e and e': 

(r,e;m\A\s,e';m) = A™ s 5 e , e , . (25) 

Substituting the operator decomposition of Eq. (|2"Tj) into 
Eq. dlTJ) and using the definition |%m(*)) = ImlVK*)) of 
Eq. (1TB)) , the time-dependent expectation value takes the 
form 

M 

A(t) = (W)\Mm)W)) + (XM(t)\A\ XM (t)) . 

m—m m i n 

(26) 
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This general expression reduces for a local operator to 

M trun 

A(t) = ( XM (t)\A\ XM (t)) + ]T £4>£(t), (27) 

m— m m j n r,s 

where 

<r(*) = E< s ' ^ ™IV>(*)> e; m) (28) 

e 

is the reduced density matrix at iteration m. 

Three comments should be made about Eqs. {27} and 
(|28[). First, these expressions are both general and ex- 
act for the real-time dynamics on the discretized chain. 
Apart from the restriction to local operators, no further 
approximations or assumptions are involved. Second, 
Eqs. (|27|) and ([28]) can be easily extended to a statisti- 
cal admixture of initial states with the statistical 
weights {wi}. This requires the simple substitutions 

(XM(t)\A\x M (t)) ^^WiixuA^lXuM) ( 29 ) 

i 

and 

m = \m)(m\ ->E^iV'i(*)><V'i(*)i (so) 

i 

in Eqs. (f2"7)) and ([2"5)> . respectively. Third, the conven- 
tional TD-NRG approach is recovered from Eqs. (|2"T|) and 
([28]) by (i) setting M = A, (ii) realizing that A x = for 
N = M, and (iii) adopting the standard NRG approxi- 
mation "H|r, e;m) « E™\r, e; m), which simplifies p™ r (t) 
to e^ E r- E T)t p f^ (m) with 

Ps?r( m ) =E<s,e;m|V'o)(V , o|r,e;m) . (31) 

e 

A natural generalization of the TD-NRG is to apply 
the NRG approximation H\r, e;m) ~ E™\r, e;m) to the 
early iterations m < M only, converting Eqs. ([27)) and 
© to 

A/ trun 

A w= E E^'Ww 

m— m m i n r,s 

+ ( XM (t)|i|XA/W> • (32) 

This equation, which constitutes one of the central results 
of this paper, interpolates between the TD-NRG, corre- 
sponding to M = A, and the exact time-dependent ex- 
pectation value, obtained for M = TO m ; n . Of course, the 
latter statement assumes an exact evaluation of \xM{t)), 
which is an impractical task for M = m m i n . As dis- 
cussed below, a proper choice of the parameter M allows 
for an improved evaluation of |xm(£)) using alternative 
methods such as the TD-DMRG or CET, with minimal 
loss of accuracy at the early iterations to which the NRG 
approximation is applied. Furthermore, by resorting to 



methods that do not rely on the special structure of the 
Wilson chain to evaluate \xM(t)), one can abandon the 
exponential decay of the hopping matrix elements be- 
yond site M, reducing thereby the discretization errors 
inherent to the Wilson chain. These principles form the 
core of the hybrid approach. We now turn to elaborate 
on the technicalities of how M is selected, the interface 
with the hybridized method, and the way in which the 
initial state jV'o) is constructed. 



C. Interface between the TD-NRG and the 
hybridized approach 

1. Hierarchy of energy scales and the time evolution of 

\XM{t)) 

To turn Eq. (|32p into an operative platform for hy- 
bridizing the TD-NRG with alternative methods of com- 
puting the real-time dynamics of |xm(^)), it is useful to 
go back to the partitioning of \tp(t)) specified in Eq. ([H|) 
and gain a deeper insight into the energy scales encoded 
into the projectors P m . Applying the operator decom- 
position of Eq. (l2~Tj) to the Hamiltonian ri, the latter is 
written as 

M 

H= E H ( TO ) + h x - (33) 

m=m m i n 

where 

H(m) = P m nP m + P m ni+ + i+nP m (34) 

and 

K = iJ f Hi+ . (35) 

It is rather easy to see that the different Hamilto- 
nian terms that appear in Eq. (|33[) generally do not 
commute with one another, i.e., [H(m),h x ] 7^ and 
^i{m),U{m')] ^ if m ml. According to the NRG 
philosophy, however, the off-diagonal terms P m HPm' 
with m ^ m! are expected to be small, as these cou- 
ple excitations on different energy scales. Consequently, 
one can approximate 'H(m) with h m = P m 'HP m to obtain 
the approximate Hamiltonian 

M 

H s» E hm + h x ■ (36) 

m=m m i n 

Evidently, Eq. (|36[) becomes exceedingly more accurate 
the smaller is M, acquiring the status of an identity for 
A I = m m ; n , since % — h x in this case. Furthermore, 
since the Hamiltonian terms h x and h m with m < M 
are confined to the subspaces projected out by and 
P m , respectively, the Hamiltonian of Eq. (|3"61 is block- 
diagonal in these subspaces with [h m , h m >] = [h m , h x ] — 
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0. This allows us to write the time-dependent state \ijj(t)) 
within this approximation as 

M 

W))= E ^ i ' hmt \^m)+e- i ' h - t \x M ) > (37) 

where \cj> m ) = P m \4>o) and \xm) = IjtflV'o) are the projec- 
tions of the initial state onto the subspaces defined by P m 
and l^f, respectively. In other terms, it suffices in this 
approximation to first project out \xm) and \4> m ) from 
the initial state, and then propagate them separately in 
time, each according to its own Hamiltonian. 

Physically, Eq. (|57|) prescribes a decomposition of the 
desired time-dependent state into independent compo- 
nents, each associated with a different time scale t m — 
1/D. m ~ A™/ 2 and evolving according to its own reduced 
Hamiltonian (either h m or h x ). In the case of spinless 
electrons the reduced Hamiltonian h m has the explicit 
form 

h m = 22E% n \l,e;rn)(l,e;m\ 

l,e 
AT-1 

+ E tniPmfl+JnPm + H.C.} , (38) 

n=m 

where /t creates an electron on the nth site of the Wilson 
chain, t n is the dimension-full hopping matrix element 
between sites n and n+1 along the chain, I runs over the 
NRG eigenstates discarded at iteration m, and -E™ de- 
notes their corresponding NRG eigenenergies. Note that 
the projection operators P m in the right-most term are 
attached in practice only to f m and f m , as all other oper- 
ators /„ with n > m do not possess any matrix element 
that takes us out of the subspace defined by P m . The 
Hamiltonian h x is nearly identical, except that the index 
m is replaced with M and the discarded states \l,e;m) 
are replaced with the NRG eigenstates retained at the 
conclusion of iteration M. In the presence of additional 
bands the Wilson orbitals /t acquire an additional flavor 
index which may label the spin cr, an orbital channel 
a, or the spin-channel tuple v = (cr, a) as in two-channel 
Kondo models (see, e.g., Ref. vm ). Other than setting 
fn fnv and adding a suitable summation over v, the 
very same equations carry over to the general multiband 
case. 



2. Physical role of the parameter M 

The Hamiltonian of Eq. can be interpreted as 

modeling a hyper-impurity with the localized configura- 
tions 1 1) and eigenenergies E™, which are tunnel-coupled 
to a chain of length N — m. The size of the impurity is 
equal to the number of states discarded at iteration m. 
Thus, the calculation of \<j> m {t)) = e~ lhmt \4>m) becomes 
exceedingly more affordable the larger is m due to the ex- 



-m/2 




Wilson chain constant hopping 



FIG. 4. A hybrid Wilson chain where the first M + 1 hopping 
matrix elements decrease according to t„ oc A"' l/2 with A > 1, 
and all further hopping matrix elements are held constant and 
equal to t_M- 



ponential reduction of the Fock space of the chain R m ,N 
attached to the hyper-impurity. We stress, however, that 
the dimension of the subspace associated with P mmin is 
comparable in size to that of the full Wilson chain, hence 
an accurate evaluation of |0m min (i)) is similar in com- 
plexity to the calculation of the full state \tf>(t)). 

There is little computational gain in implementing 
Eq. (|37|) if all components of the wave function must 
be accurately computed. Fortunately, this generally is 
not the case for the class of problems of interest, where 
|^o) is some low- lying eigenstate (typically the ground 
state) of an initial Hamiltonian rt l . Under these circum- 
stances | ?/;o) typically has only negligible overlap with the 
high-energy states of W, i.e., (0 m |0 m ) <§C 1 for the initial 
NRG iterations. Overlap becomes significant only upon 
approaching a characteristic energy scale Dm where the 
spectra of the full and the unperturbed Hamiltonians be- 
gin to notably deviate from one another. Usually this 
happens at some characteristic low-energy scale of the 
problem, e.g., the Kondo temperature Tk in case of the 
Kondo Hamiltonian. 

Consequently, the initial NRG iterations with m < M 
can be treated using further approximations such as set- 
ting t n = in Eq. (1551 . corresponding to the standard 
NRG approximation. Since the reduced density matrix 
p m r (t) requires only the matrix element (ip(t)\r, e; m), 
one can implement e *|r, e; m) instead of propagating 
\<t>m(t)} in time, which simplifies the exact result of 
Eq. (|27p to the approximate expression of Eq. (1521) . The 
computational effort can therefore be focused on evaluat- 
ing \xM(t)), which dominates the expectation value A(t). 
Most importantly, given the initial state \xm) and the ef- 
fective Hamiltonian h x generated by the NRG, \xM{t)} 
can be computed using one's method of choice. 

The physical role of the integer M, which so far served 
as a mere parameter, is now disclosed: it defines the 
NRG iteration M beyond which the time-dependent state 
should be accurately computed. Moreover, this partition- 
ing can be used to improve the discrete representation 
of the continuous bath as noted above. For example, 
one can design a hybrid chain such that all sites up to 
m = M are discretized with A > I and all further sites 
are converted to A — >• I + (see Fig. 2]). Such a chain is im- 
practical for pure NRG-based calculations, but is made 
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possible by resorting to alternative methods for tracking 
the time evolution of \xu{^))- I n this manner discretiza- 
tion errors are significantly reduced at the energy scale 
Dm, corresponding to the time scale tu — I/Dm- The 
crucial point to notice is that the reduced Hamiltonian 
h x has the effective bandwidth Dm oc A~ m / 2 -C D and 
acts on a reduced chain of length N — M . This enables 
access to long time scales of order t u ^> l/D using tech- 
niques such as the TD-DMRG or CET, which otherwise 
are restricted to far shorter times. 

The only remaining uncertainty pertains to a suitable 
choice of the iteration number M . In the absence of a 
sharp mathematical criterion, the choice of M should be 
considered on a case-by-case basis. Qualitatively, one 
expects the scale Dm to correspond to max-{To, \td\} for 
the resonant- level model [see Eq. (|72t below] , and to the 
Kondo temperature^ Tr- for the Kondo model. The case 
of an Anderson impurity is clearly more subtle, as spin 
and charge relax on different time scales^ Here a differ- 
ent optimal choice of M may apply to observables acting 
on the spin and charge sectors. 



3. Construction of the initial state |i/)o) 

So far, we have assumed the decomposition of the ini- 
tial state |-0o) according to Eq. (fl"4"|) . but did not specify 
how l^o) is obtained in practice. The construction of |?/>o) 
and its projection \xm) onto the low-energy subspace de- 
fined by 1 M depends in detail on the method hybridized 
with the TD-NRG. Our discussion below covers both the 
TD-DMRG and CET. 

We begin with the initial NRG run, which provides 
us with the low-energy Hamiltonian h l x = ^m^^m 
corresponding to the initial Hamiltonian TV . Here lf M 
denotes the projection operator onto the low-energy sub- 
space of rl % retained at the conclusion of iteration M. As 
detailed in Eq. (|55| . h % comprises of a hyper-impurity, 
a residual chain of length N — M, and a tunnel cou- 
pling between both parts of the system. In the next step 
the ground state |t/>o) of h l x is computed. In case of the 
TD-DMRG this is done using the standard DMRG al- 
gorithm^ while for the CET (for which a shorter chain 
Rm,n is employed) the Davidson method^! can be used. 
At the conclusion of this step one has the initial state 
|^o ) at hand, expressed via the kept NRG eigenstates of 

rt a 



with 



I-M- 



\ipo) =J2 c k„e\h,e;M) 



(39) 



ki,e 



Given \ipa), the state \xm) is obtained by projecting 
\ipo) onto the low-energy subspace of the full Hamiltonian 
defined by 1^-. This in turn yields 



\XM, 



= J2 b k,e\k,e;M) 



(40) 



bk.e — ^ St 

ki 



k,.k 



(41) 



where S(M) is the overlap matrix defined in Eq. (|12l) . 
This state is then propagated in time according to 
|Xm(*)) = e~ ih x*\ XM ) using either the TD-DMRG or 
CET and fed into Eq. (|32|) . As for the reduced density 
matrices p r s c ^(m) entering Eq. (j3"2")l . these are computed 
recursively from \ipo) using the standard TD-NRG algo- 
rithm^ 



III. REPEATED SWITCHINGS 

Armed with the hybrid-NRG platform, we proceed in 
this section to the second major result of this paper — the 
extension of the TD-NRG from a single quantum quench 
to repeated switching events. 

Coherent control of small nanodevices such as semi- 
conductor quantum dots or superconducting flux qubits 
can be achieved by applying gate- voltage pulses or time- 
dependent electromagnetic fields. For example, circularly 
polarized laser pulses are used to induce and control the 
spin polarization in semiconductor quantum dots4£ Al- 
ternatively, one can apply rapid changes to close-by gate 
voltages in order to control the energy levels and/or the 
tunneling rates of a quantum dot. Each of these pro- 
tocols involves repeated switchings between two distinct 
configurations of the applied fields, which we denote here- 
after by a and b. Theoretically this corresponds to peri- 
odic conversions in time between two quantum-impurity 
Hamiltonians, T-L a and 'H b , that differ in those compo- 
nents describing the isolated dot and its coupling to the 
leads. Our goal is to track the time evolution of local ex- 
pectation values in response to such a sequence of switch- 
ing events. 

To clearly formulate the problem, we assume that the 
system resides at time t = in a low-lying eigenstate of 
the Hamiltonian H a , when its Hamiltonian is abruptly 
converted from Ti a to T-i b . The system then evolves in 
time under the influence of 'H b up to time t > 0, when 
the Hamiltonian of the system is suddenly switched back 
to rl a for the duration r < t < 2r. This sequence of 
switchings is repeated periodically with a period of 2t, 
as described by the time-dependent Hamiltonian 



Hit > 0) 



H a 



(2n- 1)t < t < 2nr 
2nr < t < (2n+ l)r 



(42) 



k,c 



(here n = 0,1,2,...). Within each time interval where 
the Hamiltonian is constant, the expectation value of 
a general local operator A is given by Eqs. (|27l) and 
(1281) where, depending on the time interval in question, 
\XAf(t)) and P™ r {t) pertain either to 'H a or T-L b . Explicitly, 

|Xm(*)) is replaced with \xti i 1 )) = ^mIVK*)), where 
1q m denotes the projection operator of Eq. (JT)) written 
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with respect to NRG eigenstates of H a (a — a, b). With 
these adjustments the resulting expression is formally ex- 
act, but requires explicit knowledge of the states \tp(t)) 
and \ X $(t)). 

Consider a particular time interval (2n— l)r < t < 2nr 
in which the system evolves according to the Hamiltonian 
W a . As discussed in Sec. Ill B 21 the expectation value 
A{t) can be approximated by applying the conventional 
NRG approximation to the early iterations m < M only. 
This yields Eq. ([32]) . where |xm(*)), the eigenenergies E™ 
and E™, and the state labels s and r correspond to H a . 
The only formal modification as compared to Eq. (|32|) 
pertains to the reduced density matrix p r s °^(m), in which 
the initial state l^o) must be replaced with \ip(t)} at time 
tin-i = (2ti — 1)t. The same set of rules carry over to 
any given time interval 2nr < t < (2n + l)r, except 
that 'H a is replaced with H b , and the initial state \tpo) 
is replaced with \ip(t)) at time t^n — 2nr. This leaves 
us with those instances in time where the Hamiltonian is 
abruptly converted from T-L a to W b or vice versa. 

For concreteness let us focus on the time instance t^n = 
2n,T, when the Hamiltonian is converted from T~L a to W b . 
As described in Eq. (|14p , the state of the system can be 
formally decomposed within the time interval (2n— l)r < 
t < 2nr into 

m))=\x { ${t)) + \W {a) {t)) (43) 

with |<^ (a) (*)) = Em<M^m(*))- Here the States 

\<pm{t)) are projected according to the NRG eigenstates 
of W 1 . The time-dependent density operator is next di- 
vided into Px (t) + where 

P ( ? ) (t) = \x { M ) (t)){xti ) (t)\ (44) 

and Sp^(t) — p(t) — Px\t)- Setting t — ► t 2n and re- 
placing po = |-0 o )(-0ol -> ^(^nJX^CWI in E q- OH), the 

reduced density matrix with respect to the eigenstates of 
V. a reads 

rf£«(m) = J>, e; m|pW(t 2n )|r, e; m> 



+*P#(m), (45) 
where we have defined 

*p#(ro) = ^( S ,e;m|^ a )(i 2 „)|r,e;m) . (46) 



An equivalent expression p r s,r (m) = Xs,r{ m ) + 
5p^s}r(rn) applies to the reduced density matrix with re- 
spect to the NRG eigenstates of the Hamiltonian H b . 

In the notations of Ref. |3l| [see Eqs.(31) and (32) 
therein], the first term x^{ m ) m Eq. (1431) is of the pure 
form x ++ (m), having none of the components x^ ( m )> 
X (to), and x ( m )- Thus, one has the exact conver- 
sion 

X (b) (m) = 5 f (m) x (a) (to) S(m) , (47) 

where S(m) is the overlap matrix between the NRG 
eigenstates of H a and H b defined by Eq. JT2]). By con- 
trast, the second term Sp^(m) in Eq. (|45|) involves 
all four components <5p ++ (m), <5/^ (to,), <5p '"("t), and 
6p~~(m), and as a result lacks an explicit relation^ to 
<5p( b )(m). Similar to Eq. (|47|) , we approximate 5p^ h \m) 
by 

V b) (to) = (to) 8p {a) (to) S(m) , (48) 
leading to the compact transformation rule 

p red (b) (m) = 5 f ( m ) p rcd (a) (m) (4Q) 

The complementary transformation rule for switching 
from H b to H a reads 

p rcd (a) (m) = ^ (m) p red (b) (m) (m) _ ^ 



It should be emphasized that Eq. (|48p is an ad-hoc ap- 
proximation, whose quality is difficult to assess a-priori. 
Its accuracy is necessarily controlled by the smallness of 
(Sip^lSip^) , which in turn is reduced the smaller is M. 
Below we present numerical results demonstrating this 
point. 



In addition to the reduced density matrix, Eq. (|32[) re- 
quires explicit knowledge of the projected state, either 
IXm W) or \Xm depending on the time interval in 
question. It is therefore necessary to formulate the trans- 
formation rule of \xm (t)) upon conversion of the Hamil- 
tonian. Focusing again on the time instance t = tin-, the 
state \xm (hn)) is formally given by 



M 

1x2? (*2n)) = i^lV'^n)) - E KmI^H^)) + it M \x$(hn)) , (51) 

m— m m in 

i 

where \4>m\t2n)) and |Xm (*2n)) are projected according to the operators P a ,m and 1+ M pertaining to the Hamil- 
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tonian H a . The right-most term in Eq. (|5lT) has the exact 
representation 

U,M\X ( M ) (hn)) = S^M)\xt ) (hn)), (52) 

which follows from the fact that P a ,m\XM '(*2n)) vanishes 
by construction for all m < M. As for the remaining 
terms on the right-hand side of Eq. (|5Tj) , these require ex- 
plicit knowledge of the states \(j)m\t 2n ))- Unfortunately, 

fa) 

it is unfeasible to keep track of the states \4> m (t)) , which 
forces yet another approximation. Our strategy is to omit 
these terms altogether, with the expectation that their 
combined contribution is typically small. Naturally, the 
quality of this approximation will depend on the choice 
of M and on details of % a and H b . This leaves us with 
the approximate transformation rule 

|Xm (*2n)) ~ S\M)\x~m {t2n)) , (53) 

along with its counterpart 

|x&Wi))«S(M)|x$(i 2n+1 )). (54) 



We are now in position to summarize the proposed 
algorithm for computing A{t) — (ip(t)\A\ip(t)) for a 
sequence of switching events. We begin by rewriting 
Eq. (35J) in the form 

M trim 

m— m m i„ r,s 

+{x i M ) (t)\A\ X ^(t)) , (55) 

where a equals a or b depending on the time interval, and 
A™'a°^ is the matrix representation of A at iteration m 
with respect to the NRG eigenstates of H a . The first step 

is to evaluate pT,r ^ (m; t = 0) and |xm (* = 0)) using the 
TD-NRG methodology detailed in Ref. ED- From this 
point on these quantities are evolved in time according 
to the following set of rules: 

(i) At time t = t 2n we switch from T~L a to T~L b by setting 

p™ d{h \m;t 2n ) = S\m) p rcd ^(m;t 2n )S(m) , 

\x { M\hn))=S\M)\ X $(t 2n )). (56) 

(ii) In the time interval t 2n < t < t 2n +i, the density 
matrices and wave function are propagated in time ac- 
cording to 

\XM(t)) = e- A x«- t ^x i M ) (t2n)) , (57) 

where h b x is the projected low-energy Hamiltonian corre- 
sponding to W. b . 

(hi) At time t = t 2n +i we switch back from H to H a by 



setting 

p rcd{a) (m;t 2n+1 ) = S(m) p rcd < b > (m; t 2n+1 ) S\m) , 

\xti\h n+1 ))=S(M)\x§\t 2n+1 )}. (58) 

(iv) The density matrices and wave function are propa- 
gated in the time interval t 2n+ \ < t < t 2n +2 according 
to 

p\ ad{a) {m-t) =^W(m;t 2n+ i)e i(£ ' m - B ™ )(t - t2 " +l) , 
\ XM {t)) = e- ih W-^xti\t 2n+1 )) , (59) 

where is the projected low-energy Hamiltonian corre- 
sponding to W. a . 

We stress that four distinct approximations enter 
Eqs. (|56)) - (|59)) : (i) A discrete representation of the con- 
tinuous bath; (ii) The standard NRG approximation 
UnIt", e; m) w E™\r, e; m) that is used to propagate the 
reduced density matrices within each time interval where 
H is fixed; (hi) The transformation rule of Eq. (|48|) for 
Sp(m) at each switching event; and (iv) The omission of 
those terms that originate from the |0 m )'s in Eq. ([ST]) 
in the transformation rule for \xM(t)} at each switching 
event. While the first two approximations are general 
in nature and apply, in particular, to the TD-NRG, the 
latter two approximations are specific to the present for- 
mulation of repeated switchings. The quality of those 
approximations, which become exact for M = m m i n , can 
be tested a-posteriori by comparison to cases where ex- 
act solutions are available, as done in Sec. IVI Al for the 
noninteracting resonant-level model. 



IV. CHEBYSHEV EXPANSION TECHNIQUE 

Our presentation thus far was quite general and did 
not specify a particular method to be hybridized with the 
TD-NRG. As a proof of principle, we shall demonstrate 
in Sec. ED the hybridization of the TD-NRG with the 
CET, whose principles and implementation are reviewed 
below. 

The CET— _ — offers an accurate way to calculate the 
time evolution of an initial state |^o) under the influence 
of a general stationary hnitc-dimcnsional Hamiltonian %: 

\m) = e- im |V>o>. (60) 

The main idea of the method is to construct a stable 
numerical approximation for the time-evolution operator 
e~ lUt that is independent of the initial state |V>o) an d 
whose error can be reduced to machine precision. Its 
limitation lies in the need to explicitly store certain states 
in the course of the calculation, which limits the size of 
the Hilbert space that can be handled. 

There are different ways to expand the time-evolution 
operator, the most direct one being the conventional ex- 
pansion of the exponent in powers of H. One would like, 
however, to use an expansion that converges uniformly, 
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independent of the initial state |V>o)- A suitable choice 
are the Chebyshev polynomials, defined by the recursion 
relation 



T n +i(z) = 2zT n (z) - T„_i(z), 



(61) 



subject to the initial conditions Tq(z) = 1 and T\(z) = z. 
As is well known, the Chebyshev polynomials can be used 
to expand any function f(z) on the interval — 1 < z < 1. 
Explicitly, f(z) is expressed as an infinite scries 



/0) = ^2b n T n (z), 

n=0 

where the expansion coefficients b n are given by 
2-6 n>0 f 1 f(x)T n (x) 



(62) 



-dx. 



(63) 



The expansion in terms of Chebyshev polynomials can be 
equally applied to any function F(z) with an arbitrary 
support A m j n < z < A max using the transformation 



z' = 2 . Z Amin - L 



(64) 



which maps the interval A m i n < z < A max onto — 1 < 
z' < 1. In doing so one expands in practice the function 
f{z') = F(z). 

Using the rules laid above, the function e~ lz is ex- 
panded for A m i n < z < A max as 



e~ lz = J2 bnT n (z') 



n=0 



with 



b n>0 = 2i n e- i *J n ( 4^ 



(66a) 
(66b) 



Here J n {x) are the Bessel functions, (p equals (A ma x + 
A m in)/2, and AA = A max — A m i n . Accordingly, the time- 
evolution operator e~ lUt is expanded as 



where 



e - mt = Y j b n {t)T n {U'), 



bo(t)=e-^J (^\ 



b n>Q {t)=2i n e' lat J 1 



(67) 

(68a) 
(68b) 



Here -©max (E m i n ) is the maximal (minimal) eigenenergy 
of H, AE equals E max -E miD , a = {E max + £ min )/2, and 



H' is the "transformed" Hamiltonian 

u' = 2 n ~ En T -l. 



(69) 



Finally, applying Eq. ([6T|) to the initial state 1-00) onc 
obtains 



|^(i)>=£6n(*)|^) 



(70) 



where the infinite set of states \(j> n ) — T n (H')\tpo) obey 
the recursion relational 



2H'\ 



(71) 



subject to the initial condition | <^>o) = IV'o) an d l^i) = 

Several comments are in order. First, all time depen- 
dence is confined in Eq. (|70)) to the expansion coefficients 
b n (t) of Eqs. |68|) . which are independent of the initial 
state IV'o)- Second, the Chebyshev recursion relation of 
Eq. (fTTj) reveals the iterative nature of the calculations. 
Starting form the initial state IV'o)) one constructs all 
subsequent states \4> n ) using repeated applications of the 
"transformed" Hamiltonian %' . Note that in practice 
only two such states need be stored in memory at each 
given time. Third, since J n (x) ~ (ex/2n) n for large or- 
der n, the Chebyshev expansion converges quickly as n 
exceeds AEt. Finally, the Chebyshev expansion has the 
virtue that numerical errors are practically independent 
of i, allowing access to extremely long times. The main 
limitation of the approach, as commented above, stems 
from the size of the Hilbert space, since each of the states 
(65) \4> n ) must be constructed explicitly. In our applications 
of the approach (where \xm) serves as the initial state), 
this Hilbert space comprises of the kept NRG states at 
iteration M - typically of the order of 2 10 - and the re- 
maining chain Rm,n, whose dimension is d N ~ M (d = 2 
being the number of distinct configurations of a single 
spinless site). As a result, application of the CET is con- 
fined to rather short chains that cannot be used to access 
arbitrarily long time scales. 



SINGLE-QUENCH DYNAMICS ON A 
HYBRID CHAIN 



As already pointed out in the introduction, the Wilso- 
nian discretization of the continuous bath significantly 
influences the quench dynamics, independent of which 
finite-size approach — exact diagonalization, TD-NRG, 
or the TD-DMRG — is used to track the real-time dy- 
namics of the system. As illustrated in Fig.[TJ deviations 
from the exact continuum-limit result stem from current 
reflections at sites along the Wilson chain, caused by the 
exponentially decreasing hopping matrix elements that 
are used. These in turn produce an exponentially de- 
creasing transport velocity. 
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In this section, we present a preliminary discussion 
aimed at demonstrating the potential of substituting the 
standard Wilson chain with a hybrid chain of the type 
depicted in Fig. 2] We shall do so by investigating a 
single quantum quench in the noninteracting resonant- 
level model (RLM), which can be solved exactly on es- 
sentially any finite-size chain using exact diagonalization 
of the single-particle eigenmodes. The availability of an 
exact analytical solution for the real-time dynamics of the 
level occupancy in the continuum limit [see Eq. (50) of 
Ref. |3l| makes this model an ideal benchmark for testing 
the quality of different hybrid chains in reproducing the 
continuum-limit result. Thus, one can clearly separate 
deviations caused by the structure of the chain from those 
that originate from further approximations underlying 
the TD-DMRG or TD-NRG. Furthermore, the RLM sets 
the stage for more complicated interacting models, such 
as the interacting resonant-level model discussed below. 

Physically, the RLM describes the coupling of a single 
spinless level to a continuous band of width 2D via the 
single-particle hopping matrix element V: 



u = 



k 



e k c' k c k 



EdSd - 



V 

7K 



k 



{4<i + H.c.}. (72) 



Here d) creates an electron on the level, c\, creates a band 
electron with momentum k, and N k labels the number of 
distinct values of k. The relevant energy scales in the 
problem include the level energy Ed, along with the hy- 
bridization width To = irgV 2 , where g is the conduction- 
electron density of states at the Fermi level. 

We are interested in tracking the real-time dynamics 
of the level occupancy ndit) = (d'(t)d(t)) in response 
to a sudden quench from (E d ,Ti) to (E^,Tf). Figure [5] 
depicts the time evolution of nd(t) after the level energy 
has been abruptly shifted from E l d = — 2r to E' d = 2r 
while keeping the hybridization width fixed at Tq. We 
used the same model parameters as in Fig. [TJ but varied 
the chain structure by considering different values of the 
partitioning parameter M and different chain lengths N. 
Explicitly, as illustrated in Fig. [4j these chains comprise 
of an initial Wilson chain with exponentially decreasing 
hopping matrix elements t m cx A~ m / 2 up to m = M, 
converting to a constant hopping matrix element t m — 
t M for M <m < N. 

As can be seen in Fig. Eta), there is a systematic im- 
provement in the agreement with the continuum-limit re- 
sult upon decreasing M. Specifically, for M = 12 the de- 
viations remain quite small up to t T ~ 8, at which point 
there is a revival of ndif) which nearly reaches its original 
value nd(t — 0) for t ■ T « 9.5. In contrast to a pure Wil- 
son chain, the current velocity becomes a constant along 
the chain sites with m > M , hence current reflections 
occur only at the end of the chain^ The time at which 
the revival of n d (t) is observed is simply given by the 
round-trip time for a charge pulse that originates from 
the impurity to reach the chain end and return. Note 
that this time is shorter for M = 12 than for M = 14. 
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FIG. 5. (Color online) Comparison of the exact continuum- 
limit solution for rid{i) in the RLM (solid black line) and its 
exact numerical evaluation for various hybrid chains of the 
type depicted in Fig. [4] The full red line depicts the case 
of a pure Wilson chain of length N = 40. Panel (a): Keep- 
ing the total chain length fixed at N = 40 and varying the 
partitioning parameter M. Panel (b): For M = 12 and two 
different chain lengths N = 40 and 112. Model parameters: 



L/ = To, m 



2r , T /D = 10~ 2 , and A = 1. 



Since charge is a globally conserved quantity, true ther- 
malization and relaxation can occur only in the thermo- 
dynamic limit N — y oo. The deviation of the total charge 
from its thermalized value following the quench is simply 
given by the difference in the equilibrated charges before 
and after the quench. Since the change in total charge is 
0(1) for such a local quench, there is an <D(1/N) contri- 
bution to each reservoir degree of freedom which can be 
safely neglected in the thermodynamic limit. 

For any finite-size chain, however, this effect remains 
finite and relevant for the long-time limit. In particular, 
the continuity equation leads to charge reflections at the 
end of the tight-binding chai n 22 ' 37 such that the round- 
trip time is controlled by the chain length. Indeed, in 
Fig. [Sfb) we compare n d (t) for two hybrid chains, each 
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partitioned at M = 12. The two chains share the same 
characteristic energy scale Dm for their partitioning, but 
one has a total length of N = 40 and the other N = 112. 
While the short-time dynamics is nearly identical and 
agrees well with the continuum-limit result, the revival 
time for the long hybrid chain is pushed way beyond the 
time interval presented in Fig. [5J Thus, discretization 
errors have been nearly eliminated on the time scale of 
interest. 

This latter example vividly illustrates the potential of 
hybridizing the TD-NRG with the TD-DMRG, as the 
tight-binding chain length involved is of typical DMRG 
size. Since in our example the effective bandwidth Dm is 
of order To, the TD-NRG generates an effective Hamilto- 
nian h x whose effective bandwidth Dm can be of several 
orders of magnitude smaller than the bare conduction- 
electron bandwidth D. This in turn allows access to 
long time scales of order 1 /To 3> 1 /D which lie beyond 
the reach of pure TD-DMRG, while significantly reducing 
discretization errors that are inherent to the TD-NRG. 

A concluding word is in order about the optimal 
choices of M and the structure of the hybrid chain. The 
answers to these questions are quite difficult and can be 
expected to depend both on the model and observable of 
interest. For instance, spin and charge excitations gen- 
erally propagate with different velocities. It is therefore 
quite feasible that the optimal choices of M and N will 
dependent on the observable in question. A more system- 
atic study of the optimal hybrid chain is left for future 
research. 



VI. PERIODIC SWITCHING 
A. Periodic switching in the resonant-level model 

So far, we have stressed the potential of using hybrid 
chains by comparing their exact numerical solutions with 
the continuum-limit result for simple quenches. In this 
section, we extend our focus to the hybrid-NRG approach 
as formulated in Sec. lIIII for periodic switching. To bench- 
mark our switching algorithm, we compare it with exact 
diagonalization solutions on a pure Wilson chain, post- 
poning further discussion of the usage of hybrid chains 
and the comparison to exact continuum-limit results. 
Therefore, all calculations are performed on a standard 
Wilson chain with identical parameters, hybridizing the 
CET with the TD-NRG at a given site M. In order to 
illustrate how the effective low-energy Hamiltonian gen- 
erated by the NRG can be fed into the CET (or into any 
other method of choice for that matter) , we consider an 
extreme wide-band limit by setting the bare bandwidth 
of the RLM to D = 10 5 r . 

Using the notations introduced in Sec. IIII1 we start 
at time t = with a system that resides in the ground 
state of H a , and switch repeatedly between Ti a and H b 
after each additional time interval r. As our basic energy 
scale we set Tq = 1Q~ 5 D, whose associated time scale 




FIG. 6. (Color online) Comparison of the hybrid NRG-CET 
approach and exact diagonalization for the RLM on a Wilson 
chain with A = 1.4 and N = 84. The TD-NRG resumma- 
tion is applied up to iteration M, beyond which the CET 
is used to track the time evolution of \XM(t)). The num- 
ber of states retained in the course of the NRG is equal to 
N 3 — 1024. The resulting occupancy rid{t) is averaged over 
N z = 8 equally distributed values of the twist parameter z, 
both for the hybrid-NRG and in exact diagonalization. Model 
parameters: r a = r b = To, E b d = — — To, and r = 5/To, 
with Yq/D = 1(T 5 . Upon decreasing M, the hybrid NRG- 
CET approach gradually converges onto the exact result. 

1 /To = 10 5 /D lies many orders of magnitude beyond the 
reach of the CET when applied directly to the full Wilson 
chain. We work with a chain of hxed length N = 84 and 
the logarithmic discretization parameter A = 1.4, such 
that D N ~ r /10. 

To initiate the calculation, we first perform an NRG 
run for the initial Hamiltonian T-L a and select the ground 
state of the JVth iteration as our initial stated ji/'o) (f° r 
even N the ground state is unique). We then perform a 
second NRG run for the other Hamiltonian T-L h , during 
which we construct the overlap matrices S(m). Techni- 
cally these two NRG runs can be performed in parallel, 
in which case the overlap matrices S(m) are calculated at 
the end of each NRG iteration. Backtracking from iter- 
ation N to iteration m m j n , the reduced density matrices 
p rcd ( a )(m;0) are computed using the standard TD-NRG 
algorithm^ These density matrices, along with the pro- 
jected state \xm (0))j are then used as a seed for the 
hybrid-NRG algorithm detailed in Eqs. (f55]) -(|59 l) . 

1. Simple extension of the TD-NRG to periodic switching 

Before discussing the full hybrid approach, let us focus 
for a moment on a simple extension of the TD-NRG to 
periodic switching, corresponding to setting M — N in 
Eqs. dUD-dSU). This implies taking |xm (*)) = through- 
out the calculation. 
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Using r = 5/r and E\ = -E a d = T , Fig. El depicts 
a comparison of the TD-NRG with an exact diagonal- 
ization solution on a Wilson chain with N — M = 84. 
Since the decay rate is given by To, the system almost 
fully equilibrates before the next switching event takes 
place. As a result the deviations of the periodic TD-NRG 
from the exact numerical solution are surprisingly small. 
In particular, the periodic TD-NRG correctly tracks the 
general structure of the precise solution. Nevertheless, 
a systematic degradation in accuracy is observed upon 
going from one switching cycle to the next. This loss 
of accuracy can be quantified by the growing disconti- 
nuity Ara^ = \rid{ti + + ) — nd(U — + )| at successive 
switching events, caused by the approximation made to 
p led ^(m;t). 

Reduction of the switching time to r = 1 /To prevents 
the system from relaxing to a new equilibrium state. In 
this case high-energy excitations, which are cut off by ne- 
glecting the three additional terms Sp^ (m), 5p ^{m), 
and bp ijn), contribute significantly to the time evo- 
lution of pit). Since switching occurs on a shorter time 
scale energy is frequently pumped into the system, excit- 
ing it in such a way that the three neglected high-energy 
contributions to p Icd W (m; t) gain increasing importance 
with time. Indeed, the periodic TD-NRG result for n<i(t), 
depicted by the blue curve in Fig. [7th), shows a sizable 
accumulated error that grows systematically in time. 

2. Hybrid approach to periodic switching 

Although the simple periodic extension of the TD- 
NRG already captures the correct short-time dynamics, 
the externally driven nonequilibrium state increasingly 
deviates with time from the transient dynamics of a single 
quantum quench. Energy is locally added and removed 
from the system, which can be partially dissipated via 
heat current flowing between the impurity and the finite- 
size bath. 

To properly capture this physics, we next employ the 
hybrid approach to periodic switching. The key idea is 
to identify a suitable low-energy subspace that is large 
enough for the neglected contributions in our approxi- 
mation to be small. A trivial limit is given by setting 
M < rn min , for which the hybrid approach reproduces by 
construction the exact solution on the finite-size chain. 
Even though this limit has no practical value, it illus- 
trates the point that a reduction in M should systemat- 
ically improve the quality of the results. 

With this understanding, we extend our discussion of 
the hybrid approach to m m i n < M < N. At each time 
interval where % is fixed we evolve \Xm it)) in time using 
the CET, which is essentially exact on all time scales of 
interest. The CET is restricted, however, in the size of 
the Hilbert space one can treat, which limits the length 
of the residual chain Rm,n extending beyond site M . 
With our machines we can comfortably access up to 2 22 « 
4 x 10 6 basis states, above which the computational effort 
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FIG. 7. (Color online) Same as Fig. [6] but for r = 1/To and 
N z = 1 (i.e., no z averaging). Note the much stronger M 
dependence as compared to r = 5/To. 



rapidly becomes too exhaustive. Keeping N s = 2 10 = 
1024 states at the conclusion of each NRG iteration, this 
sets an upper limit of 12 addition Wilson sites beyond 
site M, i.e., M > N - 12 = 72. To reduce finite-size 
effects we average over N z equally distributed values of 
the twist parameter— z £ (0,1]. Below we present results 
for different values of M > 72 and N z > 1. 

Figures [6] and [7] show the time evolution of rid(t) in 
response to repeated switchings between E d = — To and 
E h d = r , keeping T a = Ti, = T fixed. Figure [5] depicts 
four successive switch events with the time interval r = 
5/To, while Fig.[7]shows ten successive switch events with 
the shorter time interval r = 1/Fo. For comparison, the 
exact time evolutions on the Wilson chain are depicted by 
the black lines, after averaging over the different values of 
z. These solutions are obtained by exact diagonalization 
of the single-particle eigenmodes of each Hamiltonian and 
using them to propagate d'(t)d(t) in time. 

Several points are noteworthy. Upon decreasing M, 
the hybrid NRG-CET curves gradually converge onto the 
exact result for the Wilson chain, both for r = 1 /Tq and 
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FIG. 8. (Color online) The time-dependent response to re- 
peated switchings between T a = To/2 and r& = To for fixed 
E2 — E h d = To/2. All other parameters are the same as in 
Fig. E] 



the hybrid-NRG and in the framework of exact diagonal- 
ization. Similar to the case of a single quantum quench, 
discretization errors become more and more pronounced 
the larger are the deviations between T-L a and H and the 
longer the time that has elapsed. Clearly, the nature of 
the perturbation (e.g., E% ^ E b d vs T a ^ r&) is also of 
importance. It should be emphasized, however, that z- 
averaging alone is insufficient for removing low-frequency 
oscillations that appear when the deviations between T-L a 
and H b are large. 

To summarize this subsection, we have demonstrated 
how the NRG can be used to systematically construct ef- 
fective low-energy Hamiltonians /i° and h°. whose band- 
width is smaller by orders of magnitude than the bare 
conduction-electron bandwidth D of the original model. 
As a result methods such as the TD-DMRG, whose ac- 
curacy is usually confined to times of order 10 2 /D, can 
be used within the hybrid platform to access exponen- 
tially long time scales. This follows from the fact that 
the new effective bandwidth w Dm oc DA~ m / 2 can be 
exponentially smaller than D. 



t = 5/r . Specifically, by M = 72 the hybrid NRG- 
CET approach essentially reproduces the exact curves. 
The largest deviations are usually found immediately fol- 
lowing a switching event, when a discontinuity generally 
develops in the occupancy produced by the hybrid NRG- 
CET approach. As noted above, the discontinuity stems 
from the approximations employed in deriving Eqs. (|56l) 
(|59|) . It decreases in size with decreasing M, reflecting a 
reduction in the accumulated weight of the |0 m (i))'s in 
favor of |xm(0) m the expansion of Eq. (fT4")l . Indeed, the 
smaller the accumulated weight of the |0 OT (i))'s the more 
accurate are the approximations employed in deriving the 
hybrid-NRG. 

For M = 84, the hybrid approach corresponds to 
the simple extension of the TD-NRG to periodic switch- 
ing. As can be seen in Fig. [TJb), the discontinuities at 
t n = tit are particularly large for the shorter switching 
time t = 1 /To since more energy is pumped into the sys- 
tem. These discontinuities are greatly reduced for the 
longer time interval r = 5/r as depicted in Fig. [6] As 
discussed above, this behavior is correlated with the fact 
that the occupancy n<j,(t) for r = 5/Fo has nearly de- 
cayed in each time interval to its new equilibrium value 
corresponding to the Hamiltonian T~i a in that time seg- 
ment. In other terms, the state of the system behaves 
as if it has effectively decayed to the new ground state, 
leaving behind only weak reminiscence of the transient 
behavior. 

Figure [8] shows the complementary response to re- 
peated switchings between two different hybridization 
widths, T a = r /2 and r h = r for fixed E% = E b d = 
To/2. As can be seen, the discontinuities at t n — nr are 
barely visible in this case, yet discretization errors do give 
rise to high-frequency wiggles that are essentially absent 
in Figs. |6]and [7] The usage of ^-averaging is essential for 
reducing these discretization effects, both at the level of 



B. Interacting resonant-level model 

Having established the accuracy of the hybrid NRG- 
CET approach for the noninteracting RLM, we pro- 
ceed to apply it to an interacting problem that lacks 
an exact reference solution. Specifically, we shall use 
the hybrid NRG-CET to track the real-time dynam- 
ics of the interacting resonant-level model (IRLM) in a 
regime not accessible to other available methods. In the 
IRLM^^Ml the resonant-level model of Eq. is 
supplemented by a local contact interaction U between 
the level and the band electrons: 

H u = u(dU-^)^-y,:clc k ,:. (73) 

^ ' k k,k' 

Here : c\c k ,:= c\c k , — 8k,k>0{~£k) stands for normal or- 
dering with respect to the filled Fermi sea. Physically, 
the contact interaction U accounts for the local capaci- 
tive coupling between the localized level eft and the band 
electrons. The total Hamiltonian of the model is given 
by 

H = H RLM + Hu , (74) 
where Hrlm represents the RLM Hamiltonian of 

Eq. G2D. 

At low energies, the IRLM is equivalent to a renormal- 
ized noninteracting RLM, both describing a phase-shifted 
Fermi liquid4i^i At resonance the model is characterized 
by the renormalized tunneling rate 

T eS *D(-£) , (75) 
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FIG. 9. (Color online) Same as Fig. H for U/D = 1 and 
E b d = -E% = r cff . Here F a = T b = r = 6 x 10~ 9 D was 
adjusted so as to maintain a fixed F c ff = 10~ J D. The black 
line shows for comparison the exact time evolution on the 
Wilson chain for U = and r„ = r b = F = 10~ 5 Z>, (taken 
from Fig. [6} ■ 



where Tq = irgV 2 is the noninteracting hybridization 
width of the RLM defined in Eg. fT2"]). a equals 25 - <5 2 , 
and 

6 = (2/tt) arctan(7T£-f7/2) (76) 

is the scattering phase shift associated with U alone 
(namely, in the absence of V). Within the NRG, T e g is 
conveniently extracted from the zero-temperature charge 
susceptibility Xc = —drid/dEd, evaluated at Ed = 0. Ex- 
plicitly, we adopt the working definition r o ff = l/(nx c ). 

In Fig. IHl we show the time evolution of nd(t) in re- 
sponse to repeated switchings between E% = —F e g and 
E h d = r e ff for U/D = 1. The bare hybridization strength 
r a = Tb = Tq = 6 x 10~ 9 D was adjusted numerically 
so as to maintain a fixed T e s = 10~ 5 D. The effect of 
a finite U is two-fold. First, it renormalizes the bare 
hybridization width from Tq to r e g, which sets the ba- 
sic time scale in the problem: l/r c g. Second, there are 
pronounced oscillations that are absent for U = 0, hav- 
ing the characteristic period t osc s» 3/r c ff. To illustrate 
this point, we have borrowed from Fig. [5] the exact time 
evolution of n<j(t) on the Wilson chain for [7 = and 
r a = r 6 = r = 1(T 5 £>. The effect of U clearly goes 
beyond just a simple renormalization of the parameters 
of the noninteracting RLM, generating new interaction- 
induced oscillations. 

To further investigate this point, we have plotted in 
Fig. [TU] the real-time dynamics in response to a single 
quantum quench from E% — — r c g to E b d — T e s, for 
different values of the Coulomb repulsion U. The bare 
hybridization width F a = Tf, = To was adjusted sep- 
arately for each value of U so as to maintain a fixed 
r c ff = 10~ 5 D. All curves begin at essentially the same 




FIG. 10. (Color online) The real-time dynamics following a 
single quantum quench from E% = —T c ft to E d — T e g , for the 
IRLM with r a = Tb = To and different values of the Coulomb 
repulsion U. Here To was adjusted separately for each value 
of U so as to maintain a fixed r e g = 10~ 5 D. All remaining 
parameters are as in Fig. [5] 



initial occupancy rid — 0.75 and decay at long time to 
rid — 0.25. Hence the low-energy fixed point, which gov- 
erns the thermodynamics, is fully determined by r e ff and 
E d . 

The intermediate dynamics, on the other hand, is quite 
sensitive to U. Starting from U = and increasing U, 
damped oscillations gradually develop with a character- 
istic period that only weakly depends on U. In contrast, 
the amplitude of the oscillations and their damping rate 
strongly depend on U. The larger is U the slower does the 
envelope function decays to zero, resulting in the emer- 
gence of two distinct time scales for large U: the period 
of oscillations t osc and the characteristic damping time 
idamp- For example, while t osc ~ 3.1/F e ff for U/D = 2, 
the corresponding damping time is of order 10/r e ff. 

To understand the origin of these interaction-induced 
oscillations, it is instructive to consider the limit of a 
large Coulomb repulsion, U — > oo. In this extreme the 
level degree of freedom <fi and the zeroth Wilson shell 
/q decouple from the rest of the chain, being confined 
to a total valence of one: S d + /J/o = 1 [this valence 
is fixed by the normal-ordered form of Ha as defined in 
Eq. (75fl ]. The corresponding subspace comprises of the 
two configurations cft\0) and /J|0), where |0) denotes the 
state in which the resonant level and the zeroth Wilson 
shell are both empty. Omitting the coupling to the rest 
of the chain, the problem has been reduced in effect to a 
2x2 matrix whose eigenenergies are 




Consequently, rid(t) must display quantum beats with 
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the frequency 



n 




(78) 



To make contact with Fig. [TU1 it is necessary to express 
the period t osc — 2tt/D, in terms of r e g = 1/(ttXc)- A 
straightforward calculation gives 



n d (E d ) 



1 l 



E d 



VIEW) 



V 2 



(79) 



resulting in T c ff = 4V/ir. 
into Eq. (J78|) then yields 



tn 



Substituting V = (7r/4)r, 



err 



r c rr y/1 + (2E d /TTT cS y 



(80) 



Finally, setting E d = T c ff to match the value used in the 
curves of Fig.HUlone obtains t osc — 3.37/r e ff, in excellent 
agreement with the period observed in Fig. 1101 

Our strong-coupling analysis clearly reveals that U en- 
ters the quantum-beat frequency only via T e g. However, 
U strongly influences the amplitude of the oscillations 
and their damping rate. In contrast to the period of os- 
cillations which is well described by U — ¥ oo, the damp- 
ing time tdamp requires a finite coupling to the rest of 
the chain. For U oo the impurity and the first Wilson 
shell decouple from the rest of the chain, resulting in co- 
herent quantum beats between the two singly occupied 
eigenstates. Once U becomes finite the system can de- 
cay incoherently to the lowest of the two singly occupied 
states through a residual coupling to the rest of the chain. 
It is this decay that determines tdamp- Since the residual 
coupling should scale as 1/U for large interactions, t<jamp 
should scale as U 2 for U D. This corresponds to a 
damping rate that falls off as 1/U 2 . 



VII. DISCUSSION AND OUTLOOK 



with and without interactions, to repeated switchings of 
its energy level and tunneling amplitude to the band. 
Such a model can be used to describe a single Coulomb- 
blockade resonance in small quantum dots in regimes 
where spin degeneracy is unimportant. 

In the absence of interactions we have critically exam- 
ined our new approach by detailed comparisons to an ex- 
act evaluation of the time evolution on the Wilson chain. 
As long as the perturbations are not too large, good accu- 
racy is maintained over a fairly large number of switching 
events. Usage of the hybrid NRG-CET greatly improves 
the accuracy in cases where large deviations develop be- 
tween the exact curve and the one produced without in- 
voking the CET, see, e.g., Fig. [71 Our present implemen- 
tation of the hybrid approach is subject to the inherent 
restriction of the CET to rather small finite-size systems. 
We expect a significant boost in accuracy and flexibility 
by combining the TD-NRG with the TD-DMRG, which is 
capable of treating far larger systems. In particular, our 
approach opens the possibility to boost the TD-DMRG 
to exponentially long time scales by using the TD-NRG 
to (i) construct an effective low-energy Hamiltonian and 
(ii) account for the short-time dynamics. This portion of 
the research is left for future work. 

The potential of the new hybrid approach was next 
demonstrated by applying it to repeated switchings in 
the IRLM. Apart from specially tuned models^ this con- 
stitutes to our knowledge the first study of such driven 
dynamics in interacting quantum impurity systems. Al- 
though equivalent at sufficiently low energies to a nonin- 
teracting RLM, the driven dynamics of the IRLM shows 
clear traces of the interaction. In particular, interaction- 
induced oscillations develop for strong Coulomb repul- 
sion, which have no equivalent in the absence of inter- 
actions. We were able to explain the period of oscil- 
lations and its weak dependence on U by invoking the 
limit U —> oo. The associated damping time scales as 
C/ 2 , and can greatly exceed the natural damping time 
l/To of the noninteracting RLM. Such a long damping 
time is strictly the effect of interactions. 



In this paper, we have extended the TD-NRG in two 
ways. First, we devised a platform for combining the 
TD-NRG with complementary methods such as the TD- 
DMRG and CET for tracking the real-time dynamics of 
quantum-impurity systems. Second, we extended the 
TD-NRG from its original realm of a single quantum 
quench to more complicated forms of driven dynamics 
where repeated switchings are applied to the system. As 
a proof of principle we combined the TD-NRG with the 
CET to compute the response of a resonant level, both 
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